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Motivated by the recent numerical evidence^ of a short-range resonating valence bond state in the 
honeycomb lattice Hubbard model, we consider Schwinger boson mean field theories of possible spin 
liquid states on honeycomb lattice. From general stability considerations the possible spin liquids 
will have gapped spinons coupled to Z2 gauge field. We apply the projective symmetry group(PSG) 
method to classify possible Z2 spin liquid states within this formalism on honeycomb lattice. It is 
found that there are only two relevant Z2 states, differed by the value of gauge flux, zero or ir, in 
the elementary hexagon. The zero-flux state is a promising candidate for the observed spin liquid 
and continuous phase transition into commensurate Neel order. We also derive the critical field 
theory for this transition, which is the well-studied 0(4) invariant theory^, and has an irrelevant 
coupling between Higgs and boson fields with cubic power of spatial derivatives as required by lattice 
symmetry. This is in sharp contrast to the conventional theory^, where such transition generically 
leads to incommensurate magnetic order. In this scenario the Z2 spin liquid could be close to a 
tricritical point. Soft boson modes will exist at seven different wave vectors. This will show up 
as low frequency dynamical spin susceptibility peaks not only at the V point (the Neel order wave 
vector) but also at Brillouin zone edge center M points and twelve other points. Some simple 
properties of the 7r-flux state are studies as well. Symmetry allowed further neighbor mean field 
ansatz are derived in Appendix which can be used in future theoretical works along this direction. 
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Quantum ground state of a spin system without any 
spontaneous symmetry breaking, the so-called spin liq- 
uid, in two or higher spatial dimensions, has been a sub- 
ject of intense research since it was first proposed more 
than thirty years ago&£. These states, sometimes called 
resonating valence bond(RVB) states, generically appear 
in two varieties, the "short-range RVB state" with a gap 
to spin-carrying excitations, and the "critical spin liquid" 
with gapless spin excitations. Recently several candidate 
materials^— have emerged for spin liquids in two spa- 
tial dimensions(2D). Interestingly they all have gapless 
spin excitations. Many parent Hamiltonians have also 
been constructed for spin liquids in 2D-ii~— . However it 
remains unclear theoretically whether a simple and nat- 
ural spin Hamiltonian, e.g. the Heisenberg model, can 
have a spin liquid ground state on some 2D lattices. For 
common bipartite 2D lattices, the square and honeycomb 
lattices, quantum Monte Carlo (QMC) 15 ' 16 and other 
calculations^— have clearly shown the long-range mag- 
netic order in the ground state of the nearest-neighbor 
Heisenberg model. Therefore frustration is usually con- 
sidered as an important ingredient for stabilizing the pu- 
tative spin liquid states. 

In an exciting paper by Meng et aiX : the half-filled 
Hubbard model on honeycomb lattice Eq. (fTJ) was care- 
fully studied by quantum Monte Carlo calculations. The 
model simply consists of hopping of electrons on nearest- 
neighbor bonds < ij > and onsite repulsion between two 
spin species labeled by a =t, i, 
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onsite repulsion U > and electron hopping t, three 
different phases were observed. With small coupling 
U/t < 3.5 the system is a semi-metal with Dirac-like 
dispersion. For large coupling 4.3 < U/t the system de- 
velops long range magnetic order. In the intermediate 
coupling region 3.5 < U/t < 4.3 a very interesting state 
with both single-particle gap and spin gap appears. Var- 
ious symmetry breaking scenarios were checked in this 
state and then ruled out. It was thus concluded that this 
state is a genuine short-range RVB state. 

This is somewhat surprising considering both weak and 
strong coupling limits. Starting from the weak coupling 
limit, with the single-particle gap develops continuously 
as observed in the calculation^, it was expected that the 
spin dynamic will either inherit the gapless nature of the 
small U semi-metal phase2a, or develop certain kind of 
spontaneous symmetry breaking. 

In the strong coupling large U —> +oo limit the 
low energy Hamiltonian is the nearest-neighbor spin- 
1/2 Hciscnbcrg antiferromagnetic(AFM) model whose 
ground state has long-range colinear Neel order— and 
must have gapless spin-wave excitations as Goldstone 
modes. Indeed a magnetic order was seen in the strong 
coupling region 4.3 < U/t in the numerical simulation 1 . 
Moreover the magnetic order parameter and spin gap 
seem to both vanish continuously at the critical point 
U/t 4.3. This raises the hope to understand the ob- 
served "short-range RVB state", at least in the large U/t 
part of the parameter range, by going from the strong 
coupling side. Although the conventional wisdom^^ is 
that such continuous quantum phase transition between 
colinear magnetic order and gapped spin liquid is impos- 
sible. 

In the strong coupling regime, with single particle gap 
much larger than the spin gap(zero in magnetic ordered 
phase), it is reasonable to describe the low energy physics 
by an effective spin- 1/2 Hamiltonian, which can be de- 
rived from the Hubbard model and should be2£ (up to 
t 4 /U 3 order) 

#spin = ~ 773") S *' S i+ Jj3 Sl ' S 3 + --- 

<ij> V / «ij>> 

(2) 

where << ij » are next-nearest-neighbor bonds. As 
the "short-range RVB" region is still close to the single- 
particle gap opening transition(Mott transition) , the spin 
Hamiltonian should be much more complex than this 
leading order Heisenberg model, i.e. have strong cou- 
plings of further neighbors and/or four and even more 
spins. Solving the exact spin model will likely not be eas- 
ier than solving the original Hubbard model. In this pa- 
per we take a different approach. Using symmetry analy- 
sis we completely classify all possible stable gapped spin 
liquid states within the Schwinger boson formalism. It 
turns out that there are only two relevant states, differed 
by the gauge invariant flux, zero or 7r, in a hexagon. Some 
signatures of these two spin liquid states will be derived 
which may be checked in numerical simulations. The 
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FIG. 1: (Color online) Mean field "phase diagram" of the zero- 
flux state. Horizontal axis is the variational parameter, ra- 
tio between next-nearest-neighbor and nearest-neighbor mean 
fiefd couplings, A2/A1. Vertical axis is the average boson den- 
sity (n)inr. The dash line (h)mf = 1 indicates the boson den- 
sity of spin-1/2 system. Solid lines are phase boundaries. The 
red solid line between the zero-flux Z2 spin liquid and the Neel 
order is a continuous transition described by the field theory 
Eq. (|22p . The vertical solid black line between the two ordered 
states is a first order transition. The blue line between the Z2 
spin liquid and the incommensurate magnetic order has yet 
to be studied but is likely a continuous transition. There is 
a very small parameter range of 0.493 < A2/A1 < 0.516 (see 
also the inset) such that a spin-1/2 system will be a gapped 
Z2 spin liquid, which is a promising explanation of the ob- 
served spin liquid 1 . The variational parameter A2/A1 can in 
principle be tuned by physical parameters. For example, as 
argued in Section llVI increase of U/t will decrease A2/A1, 
which can drive a continuous magnetic ordering transition at 
the crossing point (black dot) of the dash line and the red 
solid line. 



zero-flux state turns out to be a very promsing candi- 
date for the observed short-range RVB state. We obtain 
a mean field "phase diagram" (Fig. [1} for it in terms of 
a variational parameter, which could qualitatively agree 
with the behavior of the Hubbard model close to the mag- 
netic transtion. Our symmetry analysis fixes symmetry 
allowed forms of further neighbor mean field couplings, 
which will be useful for later theoretical studies of spin 
liquids on honeycomb lattice. 

The outline of this paper is as follows. In Section HI1 we 
briefly describe the formalism of Schwinger boson mean 
field theory. In Section IlIII we apply the projective sym- 
metry group method developed in Refill to classify all Z2 
Schwinger boson states on honeycomb lattice. Details of 
the derivation are presented in Appendix [A] Two out of 
32 possible Z2 states are particularity relevant here and 
we derive the mean field ansatz up to fourth neighbors in 
Appendix [B] In Section IIVI we study some simple prop- 
erties of the two Z2 Schwinger boson states emerged from 
the PSG analysis. And we derive the continuum field the- 
ory for the transition from the zero-flux Z2 spin liquid to 
the Neel order in Appendix [Cj Conclusions and outlook 
of further developments are summarized in Section [V] 
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II. SCHWINGER BOSON MEAN FIELD 
THEORY FOR Z 2 SPIN LIQUIDS. 

A microscopic theory of spin liquid usually in- 
volves fractionalized spin-carrying particles, the spinons, 
which are strongly coupled to certain emergent gauge 
fielc£i2£r— . It is generally believed that, when the 
spinons are gapped, the system is stable only if the gauge 
field takes discrete value s 25 ' 31 (some exotic counter- 
examples exist like the doubled Chern-Simons model of 
Levin and Wen^ 3 - but will not be considered here). The 
natural candidate of such discrete gauge field for short- 
range RVB state is the Z2(Ising) gauge theory 3 ^. Thus 
throughout this paper we will assume a Z2 spin liquid 
state on the honeycomb lattice without breaking of any 
physical symmetry. 

There are several serious problems of the Z2 spin 
liquid assumption in the context of the QMC result 1 . 
First if the magnetic ordered phase is continuously con- 
nected to a Z2 spin liquid, it will usually be non-colinear 
and incommensurate 5 , unlike the observed commensu- 
rate Neel-type order. However it will be seen later in 
this paper that this expectation is not correct on honey- 
comb lattice. Also it seems that the possibility of non- 
colinear magnetic order has not been carefully checked 
in the paper by Meng el alX. Thus we believe this ar- 
gument against a Z2 spin liquid explanation may be cir- 
cumvented. The second problem is the claim made by 
Meng et alX that topological degeneracy was not ob- 
served, while a Z2 spin liquid on a torus should have 
four-fold degenerate ground states. But it was acknowl- 
edged that their numerical method might have missed 
the degenerate ground states in other topological sec- 
tors. Despite this uncertainty we believe that it is still 
meaningful to thoroughly study the possibilities of Z2 
spin liquids on honeycomb lattice. 

Another issue for the Schwinger boson formalism is 
that it is not convenient for the description of the seem- 
ingly continuous Mott transition around U/t tzs 3.5 in 
the numerical results^. We will refrain from consider- 
ing that parameter range in this paper, and strictly limit 
ourselves in the strong coupling region with large single 
particle gap. 

To continuously evolve from a magnetic ordered state 
to a Z2 spin liquid with spin gap, a natural ap- 
proach is to decompose each spin into two bosonic 
spinons, the Schwinger boson a 5 ' 30 ' 31 . The magnetic or- 
dering transition then becomes the condensation of these 
boson o 5 ' 30 ' 31 ' 35 . And a large- N Sp(iV) generalization has 
been formulated to study the problem in a controlled 
1/N expansio n 5 ' 30 ' 31 . It is also possible to get a gapped 
Z2 spin liquid from fermionic spinons 3 - 2 - but that scenario 
will not be considered in this paper. In this paper we will 
not use the Sp(iV) language, but the PSG analysis can 
be directly applied to the large- N theory. 

In the following we briefly recall the formulation of the 
Schwinger boson mean field theory. More details can be 
found in, for example, Refi^ 5 .. 



The bosonic representation of spin Si on site i is 

a,/3 

with boson operators b, spin indices a, (3 =f, i, and Pauli 
matrices <r. For this to be a faithful representation of the 
spin system a constraint on the total boson number must 
be imposed, 

Ai = £ 6 Ua = 2S (4) 

a 

where S is the size of the spin. For spin-1/2 model, 
S = 1/2, the boson density should be unity. This hard 
constraint will be relaxed in the mean field treatment so 
it is only satisfied on average under the mean field state, 

(ni)MF = « (5) 

where (-)mf means expectation value in the mean field 
theory, and the average boson density n can also be taken 
as a parameter 35 . 

Possible mean field decouplings of Heisenberg interac- 
tion Si ■ Sj can be suggested from the operator identities 

S, • S 3 = -2A\ 3 A l3 + (l/4)n,ftj 

- ~ + « - (6) 

= - (l/4)Mj + -'/',;/>•,, = B%B l3 - A%A l3 

where A y = (1/2) (6^6^ - h^^) and B ij = 
(l/2)(b| t 6- t + b\i b ji) are both SU(2) invariant. 

A mean field theory for Heisenberg AFM model will 
generally include both A and B terms 3 ^— , 

Hmf = £;.i; r i,; - BfjBij + H.c.) + - «) 

i,j i 

+ J2(A; j A ij -B* j B ij )/J ij 

(7) 

where Aij — —Aji, B^ = B^ are complex numbers called 
the mean field ansatz, and the chemical potential /ij is 
introduced to achieve the average constraint Eq. ([5]) . For 
translationally invariant states [ii = fx are uniform. And 
Aij (B^ ) on symmetry related bonds will have the same 
magnitude. Both A and B terms have been consistently 
generalized to the theory of Sp(A^) magnets and the mean 
field Hamiltonian can be regarded as a saddle point so- 
lution of the Sp(AT) action after Hubbard-Stratonovich 
transformation 39 . Here we will not use the Sp(iV) lan- 
guage and we will regard the mean field theory as a vari- 
ational approach for general spin models even beyond 
Heisenberg model. 

The mean field Hamiltonian can be diagonalized to 
solve for boson dispersions. For small boson density k 
the bosons will be gapped. Increasing boson density will 
cause boson condensation at a critical boson density k c , 
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which corresponds to a magnetic ordering transition, and 
the details of the magnetic order can be derived from the 
structure of the boson condensates^. 

For the Heisenberg model, the mean field ansatz can 
be solved from the self-consistent equations, 

(^ij/MF = — -Aij/ Jij, (Bij)MF — —Bij/ Jij, (8) 

together with the average constraint Eq. Self- 
consistent equations for non-Heisenberg models can in 
principle be derived as well. 

As discussed in Ref<2£, for the emergent gauge theory 
to be Z2, it will need either both ansatz and B^, or 
only ansatz Aij but with geometric frustration. Nearest- 
neighbor ansatz A<ij> on honeycomb lattice is bipartite 
and will lead to a U(l) gauge theory. Since the spin 
Hamiltonian Eq. ([2]) have strong further neighbor cou- 
plings, it is natural to assume that next-nearest-neighbor 
A <<i j >> is nonzero, which is sufficient to "Higgs" the 
U(l) gauge field into Z2. 



III. PROJECTIVE SYMMETRY GROUP OF 
SCHWINGER BOSON MEAN FIELD THEORIES 
ON HONEYCOMB LATTICE 

The mean field theory Eq. (JT]) is not invariant under 
the local U(l) gauge transformations of the Schwinger 
bosons 



u j a 7 ^ u j a i 
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(9) 



where the phase <f>(J) can depend on site j. The ansatz 
will transform accordingly as 

A i:j ->■ e^^+^Uy, ->■ e'l-^'+^^B.j (10) 

However the physical spin state is gauge invariant if the 
constraint Eq. (j4} is implemented exactly. Thus differ- 
ent mean field ansatz may correspond to the same phys- 
ical state. Moreover the physical symmetries, e.g. the 
space group symmetry, may not be explicitly present 
in the mean field ansatz. And it is not straightforward 
to test whether a given mean field ansatz actually con- 
forms all the physical symmetries under the constraint 
Eq. ((3]). It was first noted by Wen and collaborators, in 
the studies of fermionic mean field theories of spin liq- 
uids, that the mean field theory should have a projective 
symmetr y 40 i 41 . Namely the mean field ansatz should be 
invariant under a combined physical symmetry group and 
gauge group operation, a projective symmetry group op- 
eration. The structure of the physical symmetry group 
constrains possible structures of this projective symmetry 
group, thus constrains possible spin liquid states. This 
idea was generalized to Schwinger boson states in Rcf. 2 ' 
and applied to triangular and kagome lattices. Here we 
will directly apply it to honeycomb lattice. More detailed 
discussion of the formalism can be found in Refill. 

The honeycomb lattice and its space group generators 
are illustrated in Fig. [5] Sites are labeled as (x, y, w) 




FIG. 2: The honeycomb lattice is shown on the left. 
Opcn(filled) circles indicate the two sublattices. ai,a2 are 
primitive vectors. For simplicity we assume the lattice con- 
stant a = ai = |aa| = 1. u,v denote the two sites within 
one unit cell. The hexagon on the right is the enlarged unit 
cell with schematic illustration of the space group generators, 
translations Ti and T2, six- fold rotation Ce, and reflection a. 



with integer x,y indicating the unit cell at xsl\ + ya.2, 
and w — u,v indicates the two sites in the unit cell. 

The space group of honeycomb lattice is generated 
by two translations T\ along ai, and T2 along a2, 
and a counter-clockwise six-fold rotation Cq around the 
hexagon center (l/3)(a! +a 2 ), and a reflection a around 
the horizontal axis through the same hexagon center. 
Their actions on the lattice are 

T\ : (x, y, w) (x + 1, y,w), w = u,v (Ha) 
T% : (x, y, w) (x, y + 1, w), w = u,v (lib) 

c . j (x,V,u) -> (-y + l,x + y-l,v) 
" (x,y,v) -> {-y,x + y,u) 



(x,y,u) -> (x + y,-y,v) 
(x,y,v) -> (x + y, -y,u) 



(lid) 



We associate a U(l) gauge group element, e l ^ x ^ de- 
pendent on site j, to each element X of the space group, 
and demand that the mean field ansatz be invariant un- 
der the combined PSG operation 



b ja -> e***l x <Mb XWa , a=T4 



(12) 



where X(j) is the image of site j under the action of 
X. The structure of the space group can be used for 
solving the allowed phase functions <j>x(j)- The solution 
is straightforward and listed in Appendix [A] In the end 
we have 

cj) Tl {x,y,w) =0, cj)T 2 (x,y,w) = p^x, 

, x(x + 2y-l) {p 7 +p 8 +pg)TT 

<PC 6 {X,y,W) =PlTT 1 , 

cp a {x,y,u) =pin[ 1- x\ + p 5 iry H , 

4> a {x,y,v) =pin[ \-x\+p 5 ny-\ . 

(13) 

with w = u,v labels sublattice, and five free integer pa- 
rameters pi,P5,P7,Ps,P9 = 0,1 mod 2. Therefore there 



FIG. 3: (Color online) The zero-flux ansatz. Left part shows 
the nearest-neighbor ansatz. Single arrow from i to j means 
A<ij> = —A<cji> — Ai > 0. All nearest-neighbor B<ij> 
must be zero according to Appendix [B] Blue dash rhombus 
is the unit cell of the mean field theory, containing two sites 
u,v. The large hexagon on the right is the enlarged mean field 
unit cell showing the next-nearest-neighbor bonds. Double 
arrow from i to j means A << ij >> = — A << ji >> = A2. All 
next-nearest-neighbor B << ij >> — +B2 are real according to 
Eq. (|B5^-Eq. pgfj l. 



FIG. 4: (Color online) The 7r-flux ansatz. Left part shows 
the nearest-neighbor ansatz. Single arrow from i to j means 
A<ij> = — A < ji > — Ai > 0. All nearest-neighbor B<i 3 > 
must be zero according to Appendix [B] Blue dash rhom- 
bus is the doubled unit cell of the mean field theory, con- 
taining four sites u,v,p,q. The large double hexagon on the 
right is the enlarged mean field unit cell showing the next- 
nearest-neighbor bonds. Double arrow from i to j means 
A«ij» — —A << j i>> = A2. All next-nearest-neighbor 
B << ij >> = ±B2 are real, with the ± signs given in Eq. (|B5a[) - 
Eq. (|B5f|) . Red thick bonds are those different from the zero- 
flux ansatz Fig. 



are at most 32 Z2 states. Requiring nonzero nearest- 
neighbor A < ij > , which is natural for strong nearest- 
neighbor Heisenberg AFM coupling, eliminates three pa- 
rameters, P5 = pi, P7 = 1 and pg — ps- If next- nearest- 
neighbor A << ij >> is also nonzero as discussed in the 
end of Section|TTl one more paremeter can be eliminated, 
ps = 1, and we are left with only one free parameter 
pi = 0, 1. So there are only two relevant Z2 states with 



b Tl (x,y,w) = 0, cj)T 2 (x,y,w) = piirx, 
x(x + 2y — 1) 7T 
~~2 2' 



, 1 \ rViV ~ 1) i 

(Pcr{x,y,u) = pi7r[ hx + yj+7r, 



(14) 



From the solutions of PSG one can construct all 
symmetry allowed mean field ansatz. The expressions 
of Aij up to fourth neighbors and By up to next- 
nearest- neighbor are listed in Appendix [B] The nearest- 
neighbor and next-nearest-neighbor A^ are also illus- 
trated in Fig. [3] and Fig. 0] for zero- and 7r-fl.ux states 
respectively. In this paper the magnitudes of nearest- 
neighbor I ^4 < i j > I and next-nearest-neighbor |j4<<y>>| 
are denoted as A\, A2 respectively. The two states 
are more intuitively distinguished by the gauge-invariant 
flux^2 in the elementary hexagon, defined as the phase 
of Aij(-Aj k )Ake(-A} m )A mn (-Ani), where the six sites 
i,j,k,i,m,n are around a hexagon. For these two states 
this flux is p\iT so the time-reversal symmetry is also sat- 
isfied. 



IV. Z 2 SPIN LIQUIDS ON HONEYCOMB 
LATTICE 



In this Section we study, within the mean field treat- 
ment, some simple properties of the two Z2 spin liquid 
states found through the PSG analysis. For simplicity we 
will only use nearest-neighbor ^4<r/> — ±A± and next- 



nearest-neighbor bonds A << ij >> = ±^2, with A\ real 
positive. The ± signs are given in Fig. [3] and Fig. 2) 
Because the spin Hamiltonian is very complicated, we 
will not compute energetics of these states and will not 
derive/solve self-consistent equations of ansatz A\,A<i. 
Instead we will treat the ratio AijA\ as a variational 
parameter and study the "phase diagram" with respect 
to it. This parameter can in principle be tuned by, for 
example, the J2/J1 ratio in the nearest-neighbor next- 
nearest-neighbor J1-J2 Heisenberg AFM model on hon- 
eycomb lattice, which is proportional to (t/U) 2 for small 
t/U [see e.g. Eq. ©]. 

Note that the J1-J2 Heisenberg model on honey- 
comb lattice has been studied within a Schwinger bo- 
son formalism by Mattsson et al£&. However only 
the nearest-neighbor A < j J> and next-nearest-neighbor 
B«ij» were used. So that theory has U(l) gauge 
field instead of Z2 and will be unstable. More recently 
Cabra et al£2- studied a J1-J2-J3 model with J3 = J2 
using Schwinger boson mean field theory. They found a 
commensurate colinear magnetic order with large J2/J1, 
which is different from the incommensurate order ob- 
tained in the present paper with large Aij A\ in the zero- 
flux state, The small J%j J\ region of phase diagram in 
Refi^S qualitatively agrees with our small A^j A\ region 
for the zero-flux state in Fig. [TJ 
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A. The Zero-flux State 

The zero-flux Z 2 spin liquid (Fig. [3]) is a promising 
candidate for the numerically observed short-range RVB 
state. It has gapped bosonic spinons coupled to Z 2 gauge 
field. And it has a continuous transition into the Neel or- 
der even with small nonzero next-nearest-neighbor mean 
field coupling A 2 , as long as A 2 < A\/2. The continuum 
field theory close to this transition is derived following the 
method in Rcf. 4 — . The effective theory shows a nontriv- 
ial coupling of bosons to the Higgs field involving cubic 
power of spatial derivatives, which allows a direct transi- 
tion from Z 2 spin liquid to Neel order. This is in contrast 
to the conventional theory of transiton between Z2 spin 
liquid and magnetic ordered stated which will generically 
give a non-colinear incommensurate magnetic order. 

The unit cell of Fig. [3] contains two sites u, v. Fourier 
transform the bosons on each sublattice (w = u,v), 

h, , — - V* r -i(kix+k 2 y) h fic^ 

u (x,y,w)a — nrr / y c L'kiua I 

V^'unit cells k 

where = k • ai j2 , the mean field Hamiltonian Eq. 
becomes, up to a constant, 

H MF-^^_ AiPi _ A *p 2 Ml2><2 J* k 

k 

(16) 

where is a four component field ^k — 

(KuV Kvp b -Kup b -k,vi) T (superscript T means 
transpose), 1 2X 2 is 2 x 2 identity matrix, Pi i 2(k) are 
2x2 anti-hermitian matrices, 

/ n +l+e' (fc i- fc 2) +e - ifc 2 \ 

fi(k) = L ei ( t ,v e ^ 1 y (17) 

and 

P 2 (k) = i[sm(k 2 ) - sin(fci) + sin(fci - k 2 )]l 2x2 . (18) 

The mean field Hamiltonian can be diagonalized by a 
Bogoliubov transformation^. The mean field dispersion 
has two branches E±, each is doubly degenerate, 

E±(k) = yjf- A\h T 2A 1 KA 2y /hf 2 - W{ftf 

(19) 

where /1 = [3 + 2 cos(fci) + 2 cos(fc 2 ) + 2 cos(fci - fc 2 )]/4, 
j 2 = 4sin(fci/2)sin(fc 2 /2)sin[(fci-fc 2 )/2], dlA 2 is the real 
part of A 2 . An example of the dispersion is shown in 
Fig. [5] 

When the dispersion is gapped, E± > 0, the average 
boson number k = (h.)mf is 

f dhdk 2 1 / | M | \j4 \ 

Since we want the system to be stable against magnetic 
ordering, we want to maximize its capability of contain- 
ing bosons. When A± and magnitude \A 2 \ are fixed, 




r t k m r 



FIG. 5: The zero-flux mean field boson dispersion E± 
Eq. (|19[) . with A2/A1 = 1/2 and average boson density 
{n)MF ~ 1 (for spin- 1/2 model), along high symmetry di- 
rections T-K-M-Y [see Fig. Oa)]. Note the very low energy 
boson modes at T point. 




FIG. 6: (Color online) (a). Hexagon is the Brillouin zone 
of the zero-flux ansatz. Central black dot is the V point 
(fci, fo) = (0,0), where boson condensation happens when 
A 2 /j4i < 1/2. When A2/A1 increases from 1/2 to +00, 
the boson condensation momenta move along the red short 
lines, 7T < jkj < 4tt/3, from T to K{K'). The V point, 
three BZ edge center M points, six T points(filled red cir- 
cle, jkj = 7r) on T-K(K') lines, and six D points(open blue 
diamond, |kj = v / 3i"/2) on F-M lines are the would-be mag- 
netic Bragg peaks for zero-flux spin liquid with A2/A1 ~ 1/2, 
namely peaks in dynamical spin susceptibility at low fre- 
quency around spin gap. (b). Hexagon is the Brillouin zone 
of the original lattice. Dash rectangle is the reduced Brillouin 
zone for the 7r-flux ansatz. In the 7r-flux state, bosons can 
condense at the momenta indicated by the filled red small 
triangles, and produce magnetic Bragg peaks with possible 
wave vectors indicated by the open blue small hexagons. 

the above boson density will be maximized if A 2 is real. 
Therefore A 2 will be assumed as real positive hereafter 
(real negative A 2 case is related to real positive case by 
a gauge transformation). 

When A 2 < A\J2 the dispersion minimum is at the T 
point, (ki,k 2 ) — (0,0), in the Brillouin zone(BZ) [see 
Fig. Oa)]. With increasing boson density the bosons 
will finally condense at the T point. Like in the trian- 
gular and kagome cas o 27 i 35 , the structure of the conden- 
sate can be determined by solving the eigenvectors of 
Eq. (1161) with zero eigenvalues at the condensation mo- 
menta. Let (£q, k 2 ) = (0,0) in Eq. (|16l) and demand (one 
of) E± to be zero, we get l/V^il = 3/2 and two eigenvec- 
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tors (1,0,0, —1) T , (0,1,1,0) T corresponding to the zero 
eigenvalues. Therefore the condensate at this momentum 
is a linear combination of these two vectors, 

(*k=(o,o)) = 2i (1, 0, 0, -1) T + z 2 (0, 1, 1, Of (21) 

Complex numbers z±, z 2 determine the orientation of 
staggered moments, as in the case of triangular lattice^ 5 -. 
Define z = (zi, z 2 ) T , then the Schwinger bosons on sub- 
lattice u(v) becomes (b a ) = z ((b a ) = ia v z*). The 
moment on sublattice u(v) is M n = (l/2)z*crz [M„ = 
(l/2)z T \-iay)(r{i(jy)z* = -(l/2)zW = -M u ]. This is 
the Neel order. 

At A2/A1 = 1/2, the minima of dispersion jump to six 
T points on the T-K(K') lines with |k| — it [BZ corner 
K(K') point has |k| = Att/3}. Further increase A 2 /A 1 
to +00 will move the minima toward the K(K') points 
[see Fig. [B£a)]. The boson condensation in this case will 
in general lead to incommensurate magnetic order. Note 
that the A2 /A\ = +00 limit is just two copies of decou- 
pled zero-flux triangular lattice Schwinger boson mean 
field theor y 27 i 35 . 

A mean field "phase diagram" in terms of the vari- 
ational parameter A2/A1 and average boson density is 
constructed as Fig. [T] There is a very small parame- 
ter range 0.493 < A2/A1 < 0.516 where the critical bo- 
son density is greater than unity, namely the spin- 1/2 
system will remain to be a gapped spin liquid. This is 
particularly promising for explaining the numerically ob- 
served transition from short-range RVB to Neel state as 
U/t is increased. Because increasing of U/t will decrease 
J2/J1 cx (t/U) 2 , and thus decrease A 2 /Ai, the spin-1/2 
system will move to the left along the dash line in Fig. [TJ 
and cross the mean field phase boundary between the 
zero-flux Z 2 spin liquid and Neel order. 

In this scenario, the spin liquid will be very close to the 
mean field tricritical point A2/A1 = 1/2 and (n)MF ~ 
1.18. Therefore the momenta of low energy bosons are 
not only the T point, but also the six T (|k| = 7r) points in 
Fig.EKa). The dispersion for A 2 /A 1 = 1/2 and (k)mf = 1 
(spin-1/2) case is drawn along high symmetry directions 
in Fig. [5] to illustrate this point. The dynamical spin 
susceptibility at low frequency around the spin gap will 
have peaks at wave vectors connecting two(can be the 
same) boson condensation momenta, these include not 
only the T point, but also three Brillouin zone edge center 
M points, and these six T points, and six other D points 
[Fig. Ha)]. 



B. Critical Field Theory for the Transition from 
Zero-flux State to Neel Order 

Considering the spatial-temporal fluctuations of the 
would-be boson condensate z in the zero-flux state close 
to the transition into Neel order, one can derive the crit- 
ical field theory. The detailed derivation is given in Ap- 



pendix [Cl The boson part of the Lagrangian reads 

C z = J d 2 r {|D T z| 2 + c 2 \D r z\ 2 + m 2 \z\ 2 
3 

+ A 3 z*[J2(ej ■ D r f]z + c.c. (22) 
3 

+ X H $ • z T (ta y )[Y,(d 3 ■ D r f]z + c.c.} 

i=i 

where r is the imaginary time, r is the spatial coodinates, 
$ ~ A 2 is the scalar Higgs field, c.c. means complex 
conjugate of the previous term, and D is the covariant 
derivative with minimal coupling to the compact U(l) 
gauge field coming from the Schwinger boson represen- 
tation. Vectors ei = (2a 2 — ai)/3, e 2 = —(& 2 +ai)/3, 
e 3 = (2ai - a 2 )/3, di = -a x , d 2 = a 2 , d 3 = a : a 2 are 
defined for convenience. The velocity c and boson mass 
m and coupling constants A3 and Xh can in principle be 
derived from the microscopic theory. Magnetic ordering 
transition happens when the mass m vanishes. 

The transformation rules of z and $ fields under space 
group symmetry can be derived from the zero-flux (pi — 
0) PSG Eq. (HI, 

Ti, T 2 : z-> z, $ -> (23a) 
C 6 : z->-a y z*, (23b) 
cr : z^-ia v z*, (23c) 

The Higgs field $ ~ A 2 transforms trivially. The La- 
grangian Eq. (f2"2"|) is invariant under the PSG. 

Note that the form of the coupling between bosons z 
and the Higgs field <& is constrained by the PSG, namely 
the microscopic lattice symmetry. It is very different 
from the typical coupling 5 which involves only one spatial 
derivative, such coupling would violate the six-fold rota- 
tion symmetry here. Naive power counting shows that 
this coupling here, with cubic power of spatial deriva- 
tives, is irrelevant, which means the Higgs field will dy- 
namically decouple from the bosons at low energy. Con- 
sidering the anomalous dimensions will not change this 
conclusion. This is why the Z 2 state here still produces 
a commensurate Neel order upon boson condensation in 
contrast to the conventional theory 5 * where it usually be- 
comes a non-colinear incommensurate order. However 
the Higgs mechanism for reducing U(l) to Z 2 is still in- 
tact, as long as the Higgs condensate $ ~ A 2 is nonzero, 
providing stability against confinement in compact U(l) 
gauge theory in 2 + 1 dimension. It would be very in- 
teresting to see if the same critical field theory can be 
reached from the Neel order side. 

At the transition point, the low energy theory is the 
0(4) invariant critical theory for the transition between 
a spiral magnet and a gapped spin liquid^— . The scaling 
properties have been studied within large- TV expansion 2,3 
and also numerically^. For example spin-spin correlations 
will have power-law scaling at large distance 

(S(0)-S(r))~|rp (24) 
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where 77 has been numerically determined^ as 77 = 
1.373(3). This can be checked with the finite-size scaling 
results of the Hubbard model when U /t is tuned to the 
magnetic ordering transition. 



looks like, up to a constant, 

/j1 4X 4 A 1 P 1 +A 2 P 2 



-4iPi - A*P 2 



fil 4 



(25) 



C. The 7r-flux State 

Now we consider the 7r-flux state in Fig. 2] The unit 
cell for the mean field theory is doubled along a 2 direction 
and contains four sites u,v,p,q. The Brillouin zone is 
halved as shown in Fig. [6jb). However we stress here 
that the physical spin state obtained from imposing the 
constraint Eq. Q on this mean field wave function has 
the original translation symmetry of honeycomb lattice, 
and this is guaranteed by the PSG. 

The mean field Hamiltonian after Fourier transform 



where 'J'k is an eight component field, Vf^ = 

14x4 is 4 x 4 identity matrix, P 12 are 4x4 anti- 
hermitian matrices, 



Pi = 7T 



1 


1 + ei 




-e- 1 + e 2 \ 

1 




(26) 



-1 + 6 



2i sin(fci) 


1 



ei + £3 





-2i sin(fci) 




l-e 2 







-1 



ei + £3 











-2i sin(fci) 




£2 - ~ ' 


2i sin(fci 



-1 
3 



) / 



with the short-hand notations t\ — e lkl , e 2 — e lk2 , £3 = 
e i(fe^-fei) ) and /ci = k • ai, fc 2 = k - (2a 2 ). Note that k' 2 is 
twice of the k 2 in previous Subsection. 

The mean field Hamiltonian can in principle be diago- 
nalized by a Bogoliubov transformation to give the mean 
field dispersion. However with A\ and A 2 both nonzero 
this is very difficult analytically. In the following we will 
set A 2 to zero and present some results for the nearest- 
neighbor ansatz. The mean field dispersion with only 
nearest-neighbor ansatz has two branches, each is four- 
fold degenerate, 

E%> (k) = ^-A?[3/4±V7(k)] (28) 

where /(k) = [3 + cos(2fci) + cos(fc 2 ) - cos(2fci - k' 2 )}/8. 
Average boson density k = (ti)mf is 

, fdhd^if \tA \jA \ 1 (29) 

The critical boson density is achieved when = 
\/3/2, and n c = 2.14 > 1. Taken at face value it means 
this state can remain quantum disordered for spin-1/2 
and even spin-1 systems. 

The bosons will condense at four momenta in the re- 
duced Brillouin zone [see Fig. Eth)], which arc k = 
±k cl = ±(fci = 7r/6,fc 2 = — 7r/3) and k = ±k c2 = 
±(fei = — 57r/6, k' 2 = — 7r/3). The condensate at each 



momentum will be 

(*k=+(ir/6,-7r/3)) = ziVi+z 2 V 2 , (30a) 

(*k=-( W /6,- W /3)) = V>iV? + W2V$, (30b) 
(*k=+(-57r/6,-7r/3)) = Z3V3 + Z3V4,, (30c) 
(*k=-(-5*/6,-,r/3)> = W 3 V 3 *+W 4 V 4 *. (30d) 

with complex cocfficcnts #1,2,3,4, ifi, 2,3,4, and the com- 
plex vectors V\,V 2 are eigenvectors of Eq. (|25|) at k c i = 
(tt/6, — 7r/3) with eigenvalue zero, and ^3,^4 are for 
k c2 = (— 57r/6, — 7r/3). The vectors Vi, 2 ,3,4 are explicitly 
given below, 

Vx = (e-"/ 12 , 0, V2 + Vi, 0, 0, -e-' i7r / 12 y / 2 + V3, 0, -1), 
^2 = (0, e-" r/12 y / 2+^, 0, -1, e - 47r/12 , 0, y/2 + V3, 0), 
y 3 = ( e 5Wi2^/ 2 + ^ o, 1, 0, 0, - e 5 Wi2 ; 0j + y/3), 

V 4 = (0, e 5 "/ 12 , 0. ^2 + V3, e 5ra / 12 y / 2 + %/3, 0, 1, 0) 

(31) 

Note that 21,2,3,4, if 1,2,3,4 may not be independent, be- 
cause one need to make sure that the number of con- 
densed bosons on every site is the same^. 

The magnetic order is complicated but will certainly 
not be the Neel order. Because bosons have to condense 
at several different momenta otherwise the condensed 
boson density(size of the magnetic moment) would be 
non-uniform on the four sublattices. Without know- 
ing the detailed condensate structure we can still de- 
termine the possible magnetic Bragg peak wavevectors, 
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which are the differences between two boson condensa- 
tion momenta. These possible magnetic Bragg peaks are 
(k%, k?,) = ±(7r/3 + m7r, — ir/3 + nir), ±(m7r, mr) with inte- 
gers m, n and are illustrated in Fig.[5f b). These momenta 
are accessible on 6x6, 12x12 and 18 x 18 lattices used 
in the quantum Monte Carlo study^. So whether this 
7r-ffux state is realized can be tested by measuring static 
spin structure factor at these momenta in the magnetic 
ordered phase. The detailed magnetic order pattern will 
be very nontrivial like that from the triangular lattice 
7r-flux stated, but will be left for future works. 

We will not study the effect of the next-nearest- 
neighbor coupling A2 in the 7r-flux state in this paper. 
We just note here that with A2/A1 — > 00, the mean 
field ansatz Fig. |4] becomes two copies of decoupled ir- 
flux states on the triangular lattice found in Refj^l. 

It would be interesting to realize this 7r-flux state in a 
simple spin model on honeycomb lattice. However for the 
nearest-neighbor Heisenberg model general argument 42 
indicates that zero-flux state will always have lower en- 
ergy than the 7r-flux state. Ring-exchange interaction 
(for the six sites around a hexagon) may favor the 7r-flux 
state 2 ^. However the natural sign of the ring-exchange 
coupling derived from the Hubbard model will actually 
favor the zero-flux state as discussed in Ref. 27 . Thus the 
7r-flux state is not likely realized in the numerical simu- 
lation of the Hubbard modeli. 



V. CONCLUSIONS 

In hope of understanding the numerical evidence of a 
short-range RVB state found by recent quantum Monte 
Carlo simulations of honeycomb lattice Hubbard model 1 , 
and the possibly continuous quantum phase transition 
from the short-range RVB to the magnetic ordered Neel 
state, we studied the Z 2 spin liquids within the Schwinger 
boson mean field theory. Applying the projective sym- 
metry group method for Schwinger boson states^ 7 , we 
completely classified possible Z 2 Schwinger boson spin 
liquid states on honeycomb lattice. Symmetry allowed 
mean field ansatz are derived for up to fourth neighbor 
couplings, which can be used for future studies of the 
Schwinger boson mean field theory. Assuming nonzero 
nearest-neighbor and next-nearest-neighbor mean field 
couplings Ai and A2, there are only two Z 2 states on 
honeycomb lattice which do not break any lattice sym- 
metry. The two states are differentiated by the gauge 
invariant flux, zero or tt, in the elementary hexagon. 

The zero-flux state is a very promising candidate for 
the numerically observed short-range RVB state. Its crit- 
ical boson density decreases from 1.18 at A2/A1 = 1/2 
to 0.516 at A2/A1 = 0, and a continuous quantum phase 
transition to Neel order will happen in this process, em- 
ulating the behavior of the numerically studied Hubbard 
model when U/t increase from below U/t — 4.3 to +00. 
The critical field theory for the phase transition to Neel 
order is an 0(4) invariant theory Eq. (|22|) , with an ir- 



relevant coupling between Higgs field and boson fields 
involving cubic power of spatial derivatives, unlike the 
conventional form of such coupling with only one spa- 
tial derivative^. Therefore it allows for a direct transiton 
from a Z2 gapped spin liquid to a Neel order. In this 
scenario the spin liquid could have soft spin fluctuations 
at not only the ordering wave vector T point, but also at 
Brillouin zone edge center M points, and six T (|k| = tt) 
points, and six other D points [see Fig. |Bfa)]. which 
can be checked by numerically calculating the dynamical 
spin susceptibility. Also the magnetic ordering transition 
will be an 0(4) invariant theory, the (finite-size) scaling 
of correlation functions can be checked against known 
results 2- — , e.g. spin-spin correlation function behaves as 
| r |-i.373 a £ i ar g e distance r. 

The 7r-flux state has the critical mean field boson den- 
sity k c f=a 2.13 (with only nearest-neighbor mean field 
couplings) well above unity. Boson condensation in the 
7r-flux state will lead to magnetic Bragg peak at several 
wave vectors as show in Fig. &[b) , including the Neel or- 
der wave vector, which can be checked in the numerical 
simulations of the magnetic ordered phase. But for en- 
ergetic reasons it is not likely realized in the Hubbard 
model. 

There are still many remaining interesting questions 
and possible future directions in this problem. (1). The 
Z2 spin liquid on a torus will have four-fold ground state 
degeneracy which was not observed in the numerical 
simulation^. It is possible that ground states in different 
topological sector actually carry different physical quan- 
tum number, e.g. quantum number with respect to six- 
fold rotation, thus not all of them were accessed in the 
simulation. It would be useful to work out these vison 
quantum numbers which can guide the search of topolog- 
ical order in the numerical work. (2). The critical field 
theory Eq. (|22|) is derived from the spin liquid side. It 
would be very interesting to start from the Neel ordered 
side and see if the same conclusion can be reached. For 
comparison to numerics it may also be useful to com- 
pute the scaling properties of other observables. Also 
the mean field tricritical point in Fig. [TJ where bosons 
condense at T and six T points, might also be of some 
interest. (3). The continuous Mott transition is not easy 
to understand with the Schwinger boson formalism, but 
is more natural in the fermionic spinon formulation. It 
may be interesting to study the Z 2 states with fermionic 
spinons, and see if a unified picture of both continuous 
Mott transition and magnetic ordering transition can be 
achieved. (4). It may be useful to derive the effective 
spin model from the Hubbard model to high orders of 
t/U , then compute energetics of the zero- flux Z2 spin liq- 
uid state and other possible states, in order to produce a 
physical (mean field) phase diagram. (5). It may also be 
useful to have a concrete simple spin model which shows 
one of these Z2 spin liquid ground states. J\ — J2 model 
may be a good example, but unfortunately has sign prob- 
lem preventing large scale quantum Monte Carlo simula- 
tions. 
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There has been a proposal of non-magnetic insulator 
state in honeycomb Hubbard model close to the metal- 
insulator transition 4 ^. Its relation to the present study 
is however unclear yet. Also in a recent paper by Xu 
and Sachdev 45 another Z2 spin liquid state was proposed 
through a different formalism. Its relation to the Z2 spin 
liquid studied here remains to be clarified. 
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Appendix A: Algebraic Solution of the Z2 PSG on 
Honeycomb Lattice 

In this Appendix we list the detailed steps for solv- 
ing the Z2 PSGs on honeycomb lattice. The algebraic 
solutions will determine all possible symmetric Z2 states 
within the Schwinger boson formalism. 

The lattice and its space group generators are de- 
scribed in Section IIIII and illustrated in Fig. [2] All inde- 
pendent commutation relations between the space group 
generators are 

T^ 1 T 2 TiT 2 ' 1 = T^ X T^T X T 2 = 

1 ^6 J 1 J 2 °6 ~ J 2 °6Jl^6 — C 6 — 

T^aTia- 1 = T^ l oT x T^ x cj- 1 = a 1 = aC 6 aC 6 = 1. 

(Al) 

where 1 is the identity element of the space group. 

For reasons discussed in Section UU we will assume the 
invariant gauge group is Z2. The generator of IGG is 

bj a -> -b ja , a =t4> Vsite j (A2) 

For each space group element X, associate a gauge 
group element [U(l) phase] exp[i0x (j)] such that the 
mean field Hamiltonian is invariant under the combined 
PSG operation 

bj a -> cxp[i(f>x(j)}b X (j)a (A3) 

Note that these phases 4>x(j) and later equations of these 
phases should be understood with implicit modulo 2tt. 

If a gauge transformation 6j s — > e % ^'bi s is applied, 
then PSG elements transform as 27 (f>x{i) 4>x{i) + 
4>(i) — 0[A _1 (i)]. Using this gauge freedom one can al- 
ways assume (on open boundary condition) 

(j) Tl (x,y,w) =0, (f>T 2 {x = 0,y,w) = (A4) 



where w = u,v labels sublattice, (x, y) labels unit cell. 

For simplicity of notations we define two forward fi- 
nite differences Aif(x,y) = f(x + l,y) — f(x,y), and 
A 2 f(x,y) = f(x,y + l)-f(x,y). 

From T x 1 T 2 T\T 2 1 = 1, convert each space group el- 
ement to its corresponding PSG element, the identity 1 
to an unknown IGG element bi a — ¥ e Wl ^bi a , we have 

Ai<t>T 2 (x,y,w) =P!tt (A5) 

with integer p\ = 0, 1 mod 2. Later used integers 
P2, 3. 4, 5, 6. 7, 8. 9 are also Z2 integers. And equations between 
them should be understood with implicit modulo 2. So- 
lution of this equation together with Eq. (|A4[) is 

4>t 2 (x, y, w) = piirx (A6) 

From this one can already conclude that the flux in the 
elementary hexagon is pin. 

At this stage there are four remaining gauge freedoms. 
These gauge transformations do not change </>Ti , <Pt 2 u P 
to IGG elements, but can be used to simplify other PSG 
elements. 

Gauge freedom I: a global phase rotation, does not 
change any PSG elements, 

b(x,y,w)a r ^ b( X y W ^ a (A7) 

This can be used to fix one of the Aij to be real positive. 
We will fix j4( , o.u)->-(o. o,v) to be real positive. 
Gauge freedom II: 

b(x,y,w)a r ^ b{x> y, w)a (A8) 

Gauge freedom III: 

bt ^ e iv(x+v) bi \ fAQ) 

u {x.y.w)a 7 L, {x.y.w)a v^^/ 

Gauge freedom IV: staggered phase rotation, 

b(x,y,u) -> e+ ^(a;,iy,u), ^(x,y,v) ~^ e ~^^(a;,y,-u) , (A10) 

From Tf 1 C & T 1 T 2 1 C 6 1 = T^CqTxCq 1 = 1 we have 
Ax<j)c e {x,y, w) = p 1 ir{x + y) + p 2 n, (Alia) 

A 2 Cpc 6 (x,y,w) = PXTTX+P3TT. (Allb) 

Its solution is 

0c 6 (x,y,w) 

x(x + 2y — 1) 
= <t>C 6 (0, 0, w) + pi-K h p 2 nx + p 3 iry 

(A12) 

If gauge freedom II is applied, p% becomes p 3 + 1, 
therefore p% can always be assumed as zero. If 
gauge freedom III is applied, p 2 becomes p 2 + 1, and 
cj>c 6 (0,0,v) becomes </>c 6 (0>0,v) + tt, therefore p 2 can 
always be assumed as zero as well. If gauge freedom 
IV is applied, §Ca (0, 0, u) becomes ^c 6 (0, 0,u) + 4> and 
4>c 6 (0, 0, v) becomes </>c 6 (0, 0, v)—(p, therefore (j>c 6 (0, 0, u) 
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and (f>c 6 (0, 0, v) can always be assumed as the same. And 
now we have exhausted all gauge freedoms. 

From T^~ 1 aT\a^ 1 = T 2 7 1 crTiT^ 1 a^ 1 = 1 we have 

Ai4> a (x,y,w) = p 4 7T, (A13a) 
A 2 <?i C r(a;,2/, w) = Pi^y + p 5 ir. (A13b) 

Its solution is 

<j> a (x,y,w) = <j> a (0,0,w) +P!iry(y- l)/2 + p 4 irx + p 5 iry 

(A14) 

From Cf = 1 we have a constraint on (f>c e (0, 0, w), 

3[0 C6 (0, 0, u) + </» Ce (0, 0, v)} + (p! + p 2 )n = p 6 7r (A15) 

From (T — 1 we have a constraint on 4> a (0, 0, w), 

</v(0, 0, u) + <p a (0, 0, v) + Tr( Pl y 2 + Pi y) = p 7 ir (A16) 

This ensures p4 — p\ mod 2 because this equation is true 
for all y. 

From aC§oC§ = 1 we have a constraint on <pc 6 (0, 0, uu) 
and 0o-(O, 0, u>), 



20 CT (O,O,«)+2(fc 6 (O,O,u) 
2^(0,0,u) + 2</) C6 (0,0,v) =p 8 7r 



(A17) 



Therefore we have 

<^(0,0,u)-<^(0,0,v) =p 9 7r (A18) 

And the solution of 0c 6 (0, 0, u;) and ^(O, 0,w) is 

</v(0,0,u) = (P7+P9)tt/ 2 mod 2tt, (A19) 
<^(0,0,u) = (p 7 -p 9 )%/2 mod 2tt, (A20) 
0c 6 (O,O,u;) = (p 7 +j5 8 +p 9 )7r/2 mod2^,(A21) 

and pi + P6 + P7 + Ps + P9 = mod 2 thus p§ can be 
eliminated. 

Considering all these constraints, P2 = P3 = 0, p± = pi, 
and pe — p\ +p 7 +Ps +P9, we will reach the final solution 
of PSG shown in the main text Eq. (|13p with only five 
free Z 2 integer parameters Pi,P5,P7,Ps>P9- 



Appendix B: Realizations of the Z2 PSG on 
Honeycomb Lattice: Mean Field Ansatz 

In this Appendix we will use the solution of PSG to 
construct symmetry allowed mean field ansatz. We will 
list the PSG allowed ansatz up to fourth neighbors of the 
honeycomb lattice. 

The algebraic solution of PSG is very general and 
usually contains many free parameters. When realized 
by a particular kind of ansatz, e.g. nearest-neighbor 
ansatz, the number of free parameter will be greatly re- 
duced because there will be further constraints on the 
PSG. For example, if Ay is nonzero, and there is a 
non-identity space group element X such that X(i) = 
j, X(j) — i, namely the bond ij maps to its inverse ji, 



then Aji — —Aij — exp[i<px(j) + i<Px (*)]Ay , therefore 
4*x{j) + ^xi}) = I" mod 2ir. All such independent non- 
identity space group elements X, which map ij to itself 
or its inverse, need to be checked. The ansatz Ay is com- 
patible with this PSG if all such checks are passed. Then 
ansatz on all symmetry related bonds can be generated 
by applying the PSG operations. 

Nearest-neighbor ansatz A<y>: Assume 
-A(o,o,u)-)-(o,o,«) = At > 0. This bond under 
a becomes its inverse (0,0, v) — > (0,0, u), then 
<^<t(0, 0, u) + 4>a{0, 0, v) = ir, therefore P7 = 1. This 
bond under T-^Cf becomes its inverse as well, then 
Ce (0, 0, i0+2&7„(l, "I, v)+2c/>c e (l, 0, tt)+&7„(l, 0, v) = 
7T, therefore P7 + Ps + P9 = 1. Also under Cq <tC§ it 
becomes its inverse, then (j>c e (1, — 1, v) + 0c 6 (0, 0, u) + 

a (O,l J tt) + <7 (O,O 1 t;) + 0Ca(O,O|«) + ^C7e(O,O,«) = 7T, 

therefore p\ + p§ + p 7 = 1. These constaints require 
P5 = Pi, P7 = 1, Ps = P9 mod 2. 

All nearest-neighbor ansatz on the lattice are 



.4 



(x,y,u)->-(x,y,v) 



A 



(x,y,u)—>(x-\-l,y—l,v) 



A 



(x,y,u)->-(x,y-l,v) 



= +A U (Bla) 
= +(-l) Piy (-l) Pl A u (Bib) 
= +Ai. (Blc) 



Next-nearest-neighbor ansatz A«y»: Assume 2nd 
neighbor A(o,o,u)— (o,i,u) is nonzero A2. This bond under 
<jCq becomes its inverse, then 4> a (0,0,u) + </>o-(0, l,u) + 
0c 6 (l, + 0c 6 (O, 0,u) = 7T, therefore pi + p 5 + p 8 = 

1. Combined with constraints from nonzero nearest- 
neighbor ansatz, this gives P5 = Pi, P7 = Ps = P9 = 1- 
So there is only one free Z 2 integer pi . 

All next-nearest-neighbor ansatz on the lattice are 



A 



*-(x,y,u)-^-(x,y-\-l,u) 
(x,y,v)—>-(x-\-l,y,v) 



A 



(x,y-\-l,u)— >(x-\-l,y,u) 



A 



(x-\-l,y,v)— >(x-\-l,y— l,v) 



-A 2 , 

-(-l) Piy A 2l 
-(-l) PlV A 2 , 
-A 2 , 



(B2a) 
(B2b) 
(B2c) 
(B2d) 



A 



(x-\-l,y,u)—^(x,y,u) 



A 



(x+l,y — l,v)—>-(x,y,v) 



+(-l)f^(-l)^A 2 ,(B2e) 
_(_l)Pi»(-l)PiA 3 .(B2f) 



With both nearest- and next-nearest-neighbor ansatz 
nonzero, there are only one free Z 2 integer pi = 0, 1 in the 
PSG solution, so there are only two different Schwinger 
mean field theories. The ansatz are pictorially shown 
in Fig. [3] and Fig. |4j They are named as the zero- 
flux(pi = 0) and 7r-flux(pi = 1) states for their different 
gauge invariant flux in a hexagon. 

Third neighbor ansatz: Assume third neighbor 
is nonzero A3. This bond under a be- 

7T, 



A 



(1,— l,«)-(0,l,tt) 



comes its inverse, then <j) a (l,— l,v) + <^<r(0, l,u) = 
therefore p? = 1. Also under C| it becomes its in- 
verse, then (f> Ce (l,0,u) + (j>c 6 (0,0,v) + (f>c e {l,0,v) + 
<f>C 6 (0, 0, u) + <i>c e (1, -1, v) + <j>c 6 (0, 1, u) = 7T, therefore 
Pi + P7 + Ps + P9 = 1- Then A3 can be nonzero only in 
the zero-flux state (pi = 0) . 

In the zero-flux state, all third neighbor ansatz on the 
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lattice are 



A 



(x-\-l,y— 1,^)— ^{x,y-\-l,u) 



{x-\-l,yiV)—¥{x,y,u) 
(x,y,v)^-(x-\-l,y,u) 



= +A 3 , 
= +A 3 , 
= +A 3 . 



(B3a) 
(B3b) 
(B3c) 





(x+l,y+l,u) 


+ (-l) PlS A 4; 


(B4a) 


A, 

^-{x^y^v 


)^(x,y-l,u) = 


+A 4 , 


(B4b) 


(x,y,v)-> 


(x-2,y+2,u) = 


+(-l) p M 4 , 


(B4c) 


(x,y,v)^> 


(x~2,y+l,u) = 


+(-l) Pl A 4 , 


(B4d) 


4/ 


)^(x,y+2,u) = 


+A 4 , 


(B4e) 


{x,y,v)^> 


(x-\-l,y — l.u) 


+ (-l) Pl M 4 . 


(B4f) 



Fourth neighbor ansatz: Assume 4th neighbor 
A( 0! o,u)^(i,i,u) is nonzero A 4 . This bond un- 
der T 2 Cq becomes its inverse, then Cei (O,O, u) + 
00,(0, l,»)+^o 8 (l,-l, «) + 0o B (-l,l, «)+^c.(l,0, «) + 
0C 6 (O, -1, u) + 0t 2 (1, 1, u) + 0t 2 (O,O, v) = 7r, therefore 
P7 + Ps + P9 = 1 • This constraint is already required by 
nonzero nearest-neighbor ansatz. 

All fourth neighbor ansatz on the lattice are 



The PSG will also impose constraints on the B^j terms 
in Eq. ([7]). For an example we consider nearest-neighbor 
Z?<y>. Assume £?(o,o,u)->(o,o,i>) is nonzero B\. This bond 
under a becomes its inverse (0,0,i>) — » (0,0, it), then 
eKp{»[0 <r (O,O,i;)-^(O,O,«)]}fli = (-l)^Bi = -B 1 = 
£?*, therefore the argument Arg(£?*/£?i) = tt mod 2tt. 
This bond under T^Cq becomes its inverse as well, then 

0Ca(l,-l,«)-0Ce(O,O,u) + ^(l,O,«)-0 Ca (l > -l > «) + 

0c 6 (l,O,u) - 0c 6 (l,O,u) - = Aig(Bf/Bi). Also un- 
der CquCq it becomes its inverse, then 0o 6 (1, — 1, i 1 ) — 
0c 6 (O,O,u) + ^(0,1, u) - 0,(0,0, «) + 0c 6 (O,O,d) - 
0C 6 (O,O,u) = = Arg(B*/Bi). These conditions im- 
ply that B\ must be zero. 

Also consider next-nearest-neighbor B << ij >> . As- 
sume next-nearest-neighbor B(o,o,«)— (o,i ,u) is nonzero 
£> 2 . This bond under aCg becomes its inverse, then 
0c 6 (l,-M) - 0C 6 (O,O,f) + ^(0,1, u) - ff (O,O,«) = 
= Arg(£?2*/i3 2 ), therefore B 2 must be real. 



All next-nearest-neighbor B 



B 



(x,y,u)— >(x,y+l,u) 



<<ij» 



= +B 2 



are 



B(x,y,v)-t(x+l,y,v) = + (-l) PlV B 2 , 
B(x,y+l,u)->(x+l,y,u.) = + (-l) Pl1 B 2 , 
B(_x+l,y,v)-Y(x+l,y—l,v) = +-^2, 



(B5a) 
(B5b) 
(B5c) 
(B5d) 



B 



(x-\-l,y : u)—>(x,y 



u) = +(-ir y (-irB 2 ,(B5e) 



the zero-flux Schwinger boson mean field Hamiltonian 
Eq. (fll)|) close to the transition to Neel order. The nota- 
tions are slightly different from Rcf. 43 . And for simplicity 
we omit the compact U(l) gauge field in the derivation, 
which can be added in the final result by promoting the 
spatial-temporal derivatives to covariant derivatives. 

Rewrite the bosons in terms of the would-be conden- 
sate modes ip at the condensation memontum k = 0, 



(x, v ,u)a = i/w(^ai + ya 2 ), 



b, ^ = i ^ - v 

u [x,y,v)ct 



& 

where ei = (2a 2 — ai)/3 is the displacement of v site 
relative to the u site in the same unit cell. 

A gradient expansion is then performed on the real 
space terms in the mean field Hamiltonian Eq. ([7]). The 
bipartite mean field couplings become, up to cubic power 
of spatial derivatives (sum over spin indices a, /3 is im- 
plicitly assumed), 

b(x,y,u)t"(x',y',v)l "(x,y,u)iP(x',y',v)'f 

= _^j 1+ar . a + <^ + (^2 lfc(r ), 
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(C2) 



where Ar = (x'ai + y'a. 2 + ei) — (xai + ya 2 ). The non- 
bipartite mean field couplings are 

b(x,y,u)^b(x',y',u)\. b[x,y,u)\P{x' ,y' ,u)^ 

= t(rl ip ua [Ar ■ d r H h \ipu/3, 



(C3) 



and 



^(x,y,v)^(x',y',v)i h(x,y,v)iP(x',y',v)t 

v ,,* r a r » i ( Ar '^) 2 , (Ar-3 r ) 3 1f , (C4) 

where Ar = (x'ai + y'a 2 ) — (xai + ya 2 ). 

Plug these relations into Eq. and use the zero-flux 
ansatz Fig. [3] with nearest-neighbor and next-nearest- 
neighbor couplings A\ > and A 2 . After collecting terms 
up to cubic power of spatial derivatives, the continuum 
limit Lagrangian C becomes 



C 



d 2 r 



(x+i,v-i,.)-f(x, v ,») = +(-l) Pl!, (-irS 2 .(B5f) + 



\/3a 2 /2 

va ) 

3 , HU^rdrf v : ; ,Uy<U 

- H 

2 4 12 



A\i> uc 



c.c. 



Appendix C: Derivation of the Continuum Field 
Theory for the Transition from Zero-flux Z2 Spin 
Liquid to Neel Order 

In this Appendix we follow the prescription of 
Sachdev 43 to derive the continuum field theory from 



A 2 (l/6)iCT^Vw[VVdj • d r ) 3 ]ip u p + c.c. 



+ a 2 (i/6)z^v: q [E(^ • 5 r) 3 ]^ + C.C.} 



i=i 



(C5) 
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where c.c. means complex conjugate of the previous term, 
V3a 2 /2 is the area of honeycomb unit cell, a = |ai| is the 
lattice constant; ei 2 ,3 are the three vectors connecting a 
u site to its nearest-neighbor v sites, 

ei = (2a 2 -ai)/3, e 2 = -(a 2 +ai)/3, e 3 = (2ai-a 2 )/3, 

(C6) 

and we also define for convenience 



di = ai, d 2 = a 2 , d 3 = a : a 2 . 



(C7) 



Note that many terms are canceled due to the geometry, 
especially the first derivative terms from the Ai term 
cancel because X^=i(dj ■ d r ) — 0. 

Define two fields z and II from linear combinations of 

ip u and ip v , 

Z a = (ipua + 4>va)/2, II Q = (lp ua - 1p va )/2 (C8) 

Plug this into Eq. (|C5I) . the Lagrangian becomes (spin 
indices a, (3 are omitted), 



J V3a 2 /2l dr dr 
+ (2/i - 2,A l )z*z + (2/i + 3Ax)n*n 
+ a 2 (A 1 /3)d r z* -drz + c.c. 

3 

+ {A l /12)z*^{e j -d T f]z + c.c. 

3=1 

3 

+ (A 2 li)z T {ia y )[Y J {d ■ d T f]z + c.c. 

3 = 1 



(C9) 



Note that terms involving both field II and spatial deriva- 



tives are omitted, as they will generate terms in the ef- 
fective Lagrangian of z with fourth or higher power of 
spatial derivatives, and the following identity has been 
used, 



5](e l .a r ) 2 = (2/3)a 2 9 1 2 



(CIO) 



i=i 



Integrate out the field II with large gap 2/j, + 3Ai , we 
get the effective Lagrangian for z 



C z = d 2 r 



(2/i + 3Ai) v / 3a 2 



d T z* ■ d T z 



+ 2 -^d r z*. dr z + 2 ^- 3A ^z*z 
3V3 V3a 2 



+ 



A! 



6V3a 2 



J^- z T iia y ){ j2(d r d r r}z+c.c. 

3=1 



E(ej ■ d r f}z + c.c. 



(Cll) 



3 = 1 



The critical point is A\j fx = 2/3 consistent with the mean 
field solution. The critical boson velocity is proportional 
to A\. After a proper rescaling of t the Lagrangian can 
be cast into the simple form of Eq. (|2"2"T) . Note that A 2 
plays the role of the Higgs field. 
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